#include<../data_tools/usefull_routines.c>

int main(int argc, char *argv[])
{
  int i,npuntos;
  double R;
  double *dfdx,*d2fdx2;
  char *potentialfile,outfiles[100];
  potentialfile = argv[1];
  npuntos = contadorParticulas(potentialfile);
  printf("\ncomputing with %d points\n\n",npuntos);

  double x[npuntos],fx[npuntos],*kappa,*logx;
  double xi[npuntos],fxi[npuntos];

  logx = (double *)malloc((size_t)(npuntos)*sizeof(double));
  dfdx = (double *)malloc((size_t)(npuntos-4)*sizeof(double));
  d2fdx2 = (double *)malloc((size_t)(npuntos-4)*sizeof(double));
  kappa = (double *)malloc((size_t)(npuntos-4)*sizeof(double));

  FILE *potfile;
  potfile = fopen(potentialfile,"r");
  for(i=0; i<npuntos; i++)
    {
      fscanf(potfile,"%lf %lf",&xi[i],&fxi[i]);
      //logx[i] = log10(x[i]);
    }

  gslSmooth(xi,fxi,npuntos);

  FILE *potfileS;
  potfileS = fopen("pot_suavisado.dat","r");
  for(i=0; i<npuntos; i++)
    {
      fscanf(potfileS,"%lf %lf",&x[i],&fx[i]);
      logx[i] = log10(x[i]);
    }

  //derivateLogjuank(x,fx,npuntos);
  derivateLog(x,fx,dfdx,d2fdx2,npuntos);
  //derivate(x,fx,npuntos);

  FILE *kappafile;
  sprintf(outfiles,"kappa_%s",potentialfile);
  kappafile = fopen(outfiles,"w");
  for(i=2; i<npuntos-2; i++)
    {
      //R = logx[i];
      R = x[i];
      // printf("%lf %lf\n",x[i+1]-x[i],log10(x[i+1])-log10(x[i]));
      dfdx[i-2] = dfdx[i-2]/x[i];
      d2fdx2[i-2] = (1.0/(x[i]*x[i]))*d2fdx2[i-2]-(1.0/x[i])*dfdx[i-2];
      //kappa[i-2] = sqrt(((3.0/x[i])*dfdx[i-2]) + (d2fdx2[i-2]));
      kappa[i-2] = sqrt( (3.0/(R*R))*dfdx[i-2] + (1.0/(R*R))*d2fdx2[i-2]  );
      fprintf(kappafile,"%e %e %e %e %e %e %e\n",x[i],logx[i],fx[i],dfdx[i-2],d2fdx2[i-2],kappa[i-2],kappa[i-2]*kappa[i-2]);
    }

  free(dfdx);
  free(d2fdx2);
  fclose(kappafile);
  fclose(potfile);
  fclose(potfileS);


  return 0;
}
